library(readr)
library(rms)
library(bkmr)
library(mice)
library(forestplot)
library(tableone)
library(corrplot)
library(ggplot2)
setwd('/home/jar/zxjar/paper_helper/userdata/0/68/')
demo <- read_csv("demo.csv")
################################## table S3 ####################################
#参与的变量
myVars1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
nonvar1 <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
tab1 <- CreateTableOne(vars = myVars1, strata =c("Hyperuricemia")  , data = demo)
tabMat1 <- print(tab1, nonnormal = nonvar1,test = TRUE)
tabMat1 <- tabMat1[,c(-4)]
sink("BaselineS3.txt")
tabMat1
sink()

######################## 邻苯二甲酸酯相关性热图 图2 ############################
demo2 <- demo[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')]
data<-cor(demo2,method="spearman")
Cairo::CairoTIFF(file="fig2.tiff", width=8, height=8,units="in",dpi=150)
corrplot(data, method = "number", type = "upper",order = 'hclust',
         tl.col = "black", tl.cex = 0.8, tl.srt = 45,tl.pos = "lt")

corrplot(data, method = "circle", type = "lower",order = 'hclust',
         tl.col = "n", tl.cex = 0.8, tl.pos = "n",
         add = T)
dev.off()


################################## RCS 图S2 ####################################
demoS2 <- demo
demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')] <- 
  log(demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')])

ddist<-datadist(demoS2)
options(datadist="ddist")
RCSX <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
RCSNK <- c(3,3,3,3,3,5,3,3,3,3)
i <- 1
Cairo::CairoTIFF(file="figS2.tiff", width=8, height=8,units="in",dpi=150)
par(mfrow=c(2,5))
for (i in 1:10) {
  RCSX1 <- RCSX[c(-i)]
  
  paste0(RCSX1,collapse = '+')
  fml <- as.formula(paste0('Hyperuricemiaint ~rcs(',RCSX[c(i)],',nk=',RCSNK[c(i)],')+',paste0(RCSX1,collapse = '+'),
                           '+Age + Sex + race + edu + FPL + activity + smoke + drink + Hypertension + CKD + urinecreatinine + BMI'))
  fit <-lrm(fml,data=demoS2,x=FALSE)
  pred_OR<-Predict(fit,name=RCSX[c(i)],ref.zero=TRUE,fun=exp)
  
  ylim.bot<-min(pred_OR[,"lower"])
  ylim.top<-max(pred_OR[,"upper"])
  plot(pred_OR[,1],pred_OR[,"yhat"], 
       xlab = RCSX[c(i)],ylab = paste0("OR"),
       type = "l",ylim = c(ylim.bot,ylim.top),
       col="red",lwd=2)+
    lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)+
    lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)+
    lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) +
    points(1,1,pch=16,cex=1.2)#需要设置参考值,
}
dev.off()


################################## RCS 图S3 ####################################
demoS2 <- demo
demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')] <- 
  log(demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')])

ddist<-datadist(demoS2)
options(datadist="ddist")
RCSX <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
RCSNK <- c(3,3,3,3,3,5,3,3,3,3)
i <- 1
Cairo::CairoTIFF(file="figS3.tiff", width=8, height=8,units="in",dpi=150)
par(mfrow=c(2,5))
for (i in 1:10) {
  RCSX1 <- RCSX[c(-i)]
  
  paste0(RCSX1,collapse = '+')
  fml <- as.formula(paste0('LBDSUASI ~rcs(',RCSX[c(i)],',nk=',RCSNK[c(i)],')+',paste0(RCSX1,collapse = '+'),
                           '+Age + Sex + race + edu + FPL + activity + smoke + drink + Hypertension + CKD + urinecreatinine + BMI'))
  fit <-lrm(fml,data=demoS2,x=FALSE)
  pred_OR<-Predict(fit,name=RCSX[c(i)],ref.zero=TRUE,fun=exp)
  
  ylim.bot<-min(pred_OR[,"lower"])
  ylim.top<-max(pred_OR[,"upper"])
  plot(pred_OR[,1],pred_OR[,"yhat"], 
       xlab = RCSX[c(i)],ylab = paste0("OR"),
       type = "l",ylim = c(ylim.bot,ylim.top),
       col="red",lwd=2)+
    lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)+
    lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)+
    lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) +
    points(1,1,pch=16,cex=1.2)#需要设置参考值,
}



############################## bkmr单变量暴漏 图S4 #############################
demo <- read_csv("demo.csv")
demoS5 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS5 <- demoS5[c(1:200),]
y <- demoS5$Hyperuricemiaint
# Z <- demoS5[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
#                'MCOP','MCPP','MEP','MBzP')]

Z <- demoS5[,c('MECPP','MnBP')]
X <- demoS5[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:2){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
Cairo::CairoTIFF(file="figS4.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)
dev.off()



###################### bkmr三个分位数的暴露-应答趋势 图S5 ######################
demo <- read_csv("demo.csv")
demo <- data.frame(demo)
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 100, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
pred.resp.bivar <- PredictorResponseBivar(fit =fitkm , min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivar, 
  Z = Z, qs = c(0.1, 0.5, 0.9))


Cairo::CairoTIFF(file="figS5.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
dev.off()


##################### bkmr三个分位数的单次暴露风险比较 图S6 ####################
demo <- read_csv("demo.csv")
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
risks.singvar <- SingVarRiskSummaries(fit = fitkm, y = y, Z = Z, X = NULL, 
                                      qs.diff = c(0.25, 0.75),
                                      method = "exact")
Cairo::CairoTIFF(file="figS6.tiff", width=8, height=8,units="in",dpi=150)
ggplot(risks.singvar, aes(variable, est, ymin = est - 1.96*sd, 
                          ymax = est + 1.96*sd, col = q.fixed)) + 
  geom_pointrange(position = position_dodge(width = 0.75)) + 
  coord_flip()
dev.off()


############################# bkmr总体效应 图S7B ##############################
demo <- read_csv("demo.csv")
demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
                 'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
                 'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
                 'Hyperuricemiaint','LBDSUASI')]
demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
               'MCOP','MCPP','MEP','MBzP')]
X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
               'drink','Hypertension','CKD','urinecreatinine','BMI')]
for(i in 1:10){
  Z[,i] <- log(Z[,i])
}
set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
risks.overall = OverallRiskSummaries(fit=fitkm,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="figS7B.tiff", width=8, height=8,units="in",dpi=150)
ggplot(risks.overall,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='yellow')+
  geom_pointrange()
dev.off()


################################# 森林图 图3 ###################################
demo3 <- demo
demo3$MECPP_quantile <- cut(demo3$MECPP,breaks = quantile(demo3$MECPP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))
demo3$MnBP_quantile <- cut(demo3$MnBP,breaks = quantile(demo3$MnBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEHHP_quantile <- cut(demo3$MEHHP,breaks = quantile(demo3$MEHHP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEOHP_quantile <- cut(demo3$MEOHP,breaks = quantile(demo3$MEOHP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MiBP_quantile <- cut(demo3$MiBP,breaks = quantile(demo3$MiBP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$cxMiNP_quantile <- cut(demo3$MiNP,breaks = quantile(demo3$MiNP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MCOP_quantile <- cut(demo3$MCOP,breaks = quantile(demo3$MCOP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MCPP_quantile <- cut(demo3$MCPP,breaks = quantile(demo3$MCPP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MEP_quantile <- cut(demo3$MEP,breaks = quantile(demo3$MEP),labels = c('Q1', 'Q2', 'Q3', 'Q4'))  
demo3$MBzP_quantile <- cut(demo3$MBzP,breaks = quantile(demo3$MBzP),labels = c('Q1', 'Q2', 'Q3', 'Q4')) 
varis<- c('MECPP_quantile','MnBP_quantile','MEHHP_quantile',
          'MEOHP_quantile','MiBP_quantile','cxMiNP_quantile','MCOP_quantile',
          'MCPP_quantile','MEP_quantile','MBzP_quantile')
sinLogistic <- function(x){
  my_list<-list(c(1))
  for (variable in x) {
    len <- length(my_list)
    f<- as.formula(paste0('Hyperuricemiaint ~ ',variable,'+urinecreatinine'))
    glmMECPP<-glm(f,data = demo3,family =  binomial )
    my_list[[len+1]] <- glmMECPP
  }
  return(my_list)
}
allSin <-sinLogistic(varis)

my_listF <- list(c(1))
for (variable in c(2:11)) {
  print(variable)
  len <- length(my_listF)
  fit.result<-summary(allSin[[variable]])
  df1<-fit.result$coefficients
  df2<-confint(allSin[[variable]])
  df3<-cbind(df1,df2)
  df4<-data.frame(df3[-1,c(1,4,5,6)])
  df4$Var<-rownames(df4)
  colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
  df5<-df4[,c(5,1,2,3,4)]
  df5$OR_mean<-df5$OR
  df5$OR<-paste0(round(df5$OR,2),
                 "(",
                 round(df5$OR_1,2),
                 "~",
                 round(df5$OR_2,2),
                 ")")
  df5$Pvalue<-round(df5$Pvalue,3)
  my_listF[[len+1]] <- df5
}
allfros <- dplyr::bind_rows(list(my_listF[[2]],my_listF[[3]],my_listF[[4]],my_listF[[5]],my_listF[[6]],my_listF[[7]],my_listF[[8]],my_listF[[9]],my_listF[[10]],my_listF[[11]]))
fp<-allfros
Cairo::CairoTIFF(file="fig3.tiff", width=8, height=8,units="in",dpi=150)
forestplot(labeltext=as.matrix(fp[,1:3]),
           mean=fp$OR_mean,
           lower=fp$OR_1,
           upper=fp$OR_2,
           zero=0.25,
           boxsize=0.2,
           graph.pos=2)

dev.off()










